#include "mex.h"
#include <math.h>

double getgaussian()
{
    double result=0;
    double x, y;
    double tmp, randmax=(double) RAND_MAX;
    
    do
    {
        x = 2.0*rand()/randmax-1;
        y = 2.0*rand()/randmax-1;
        tmp = x*x+y*y;
    }
    while(tmp>=1.0);
    
    result = x*sqrt(-2*log(tmp)/tmp);

    return result;
}


double summation(double ptr[],int term)
{
    int i;
    double sum=0;
    for(i=0; i<term; i++)
        sum += ptr[i];
    
    return sum; 
}

static void mcforpart2(double p[], double n1[], double n2[], int r, int c)
{
        int i,j,k;
        double dt=0.5/5000;

        for (j=0; j<c; j++)
		{
            p[j*5001]=0;
			for (i=1; i<5001; i++)
            {
                p[j*5001+i]=p[j*5001+i-1]+n1[j*5001+i-1]*n2[j*5000+i-1]*sqrt(dt);
            }
		}


}



//  static void variation(long Npaths)
//  {
//      int i;
//      for (i=0,i<Npaths,i++)
//
//  }


void mexFunction(int nlhs,mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
        double *p, *n1, *n2;
        int row,column; 
        
        row=mxGetM(prhs[0]); //OU process inside, 5001 rows 
        column=mxGetN(prhs[0]);

        plhs[0]=mxCreateDoubleMatrix(5001,column,mxREAL);
        p=mxGetPr(plhs[0]);
        n1=mxGetPr(prhs[0]); 
        n2=mxGetPr(prhs[1]); //W inside, 5000 rows

        mcforpart2(p,n1,n2,row,column);


}
